//
//  MPIMatrixOperations.h
//  PIDAMRX
//
//  Created by Patrick Zulian on 4/19/11.
//  Copyright 2011 ICS - USI - Lugano. All rights reserved.
//

#ifndef PIDAMRX_MPIMatrixOperations_H
#define PIDAMRX_MPIMatrixOperations_H

#include <mpi.h>
#include "Communicator.h"
#include "CartesianTopology.h"


namespace nsesolver {
    namespace mpi {
        namespace la {
            
            
            template< class Matrix >
            class PLAlgebra {
            private:
                CartesianTopology cart;
            public:
                typedef typename nsesolver_traits< Matrix >::Scalar Scalar;
                
                explicit
                PLAlgebra(const CartesianTopology &cart)
                : cart(cart)
                {
                    
                }
                
                
                
                /**
                 * \param m1 is a distributed matrix. 
                 * \param m2 is required to be replicated for each row block in y direction
                 *  Example:
                 *  m1                m2 
                 *  _____________    _______
                 *  |   |   |   |
                 *  | P1| P3| P5|    P1 & P2 
                 *  |   |   |   |    -------
                 *  -------------    P3 & P4
                 *  | P2| P4| P6|    -------
                 *  |   |   |   |    P5 & P6       
                 *  -------------    -------
                 *
                 * \return the result of the multiplication
                 */
                const Matrix mul(Matrix &m1, const Matrix &m2 )
                {
                    
                    Communicator &comm(cart.getCommunicator());
                    Matrix result(m2.rows(), m2.columns());
                    
                    if (cart.getNumberOfDims() == 1 && m2.columns() == 1) {
                        //only matrix vector mul for now A * v
                        for (int r = 0; r < m2.rows(); r++) {
                            result[r] = m1.row(r) * m2;
                        }
                        comm.allReduce(&result[0], result.size(), MPI_SUM);
                    } else {
                        NSE_SOLVER_ASSERT( false ); //TODO
                    }
                    
                    
                    
                    return result;
                    
                } 
                
                inline
                const Scalar normInf(const Matrix &m)
                {
                    using std::max;
                    Scalar result = 0;
                    Communicator &comm(cart.getCommunicator());
                    
                    
                    for (int r = 0; r < m.rows(); r++) {
                        result = max( m.row(r).sumAbs(), result );
                    }
                    
                    comm.allReduce(&result, 1, MPI_MAX);
                    return result;
                }
                
                inline
                const Scalar norm2(const Matrix &m1)
                {
                    
                    NSE_SOLVER_ASSERT(m1.rows() == 1 || m1.columns() == 1);
                    Communicator &comm(cart.getCommunicator());
                    
                    
                    Scalar result = 0;
                    
                    for (int i = 0; i < m1.size(); i++) {
                        result += m1[i]*m1[i]; 
                    }
                    
                    comm.allReduce(&result, 1, MPI_SUM);
                    return result;
                }
                
                inline
                const Scalar distance(const Matrix &m1, const Matrix &m2)
                {
                    return norm2(minus(m1,m2));
                }
                
                inline
                const Matrix minus(const Matrix &m1, const Matrix &m2)
                {
                    return m1 - m2;
                }
                
            };
            
        }
        
    }
}

#endif //PIDAMRX_MPIMatrixOperations_H

